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Abstract 

In this paper we present a two-dimensional, time dependent, nu- 
merical simulation of a reconnection current layer in incompressible 
resistive magnetohydrodynamics with uniform resistivity in the limit 
of very large Lundquist numbers. We use realistic boundary condi- 
tions derived consistently from the outside magnetic field, and we also 
take into account the effect of the back pressure from flow into the the 
separatrix region. We find that within a few Alfven times the system 
evolves from an arbitrary initial state to a steady state consistent with 
the Sweet-Parker model, even if the initial state is Petschek-like. 

PACS Numbers: 52.30. Jb, 96.60.Rd, 47.15.Cb. 



1 Introduction 

Magnetic reconnection is of great interest in many space and laboratory 
plasmas [1], |2|] , and has been studied extensively for more than four decades. 
The most important question is that of the reconnection rate. The pro- 
cess of magnetic reconnection, is so complex, however, that this question 
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is still not completely resolved, even within the simplest possible canoni- 
cal model: two-dimensional (2D) incompressible resistive magnet ohydro dy- 
namics (MHD) with uniform resistivity rj in the limit of S — > oo (where 
S = VaL/t] is the global Lundquist number, L being the half-length of the 
reconnection layer). Historically, there were two drastically different esti- 
mates for the reconnection rate: the Sweet-Parker model [|3], ^J] gave a rather 
slow reconnection rate (E SP ~ S^ 1 / 2 ), while the Petschek || model gave any 
reconnection rate in the range from E$p up to the fast maximum Petschek 
rate Epetschek ~ 1/ log S. Up until the present it was still unclear whether 
Petschek-like reconnection faster than Sweet-Parker reconnection is possible. 
Biskamp's simulations |J are very persuasive that, in resistive MHD, the rate 
is generally that of Sweet-Parker. Still, his simulations are for S in the range 
of a few thousand, and his boundary conditions are somewhat tailored to the 
reconnection rate he desires, the strength of the field and the length of layer 
adjusting to yield the Sweet-Parker rate. Thus, a more systematic boundary 
layer analysis is desirable to really settle the question. In particular, one 
needs an elaborate and detailed picture of the reconnection current layer - 
namely, a picture that features a realistic model for the variation of the out- 
side magnetic field along the layer, and realistic 2D profiles of the plasma 
parameters inside the layer. 

The development of such a framework is the main goal of the present pa- 
per. We believe that the methods developed in this paper are rather universal 
and can be applied to a very broad class of reconnecting systems that include 
more realistic physics. However, for definiteness and clarity we keep in mind 
a particular global geometry, that presented in Fig. [I] (although we do not 
use it explicitly in our present analysis). This Figure shows the situation 
somewhere in the middle of the process of merging of two plasma cylinders. 
Regions I and II are ideal MHD regions: regions I represent unreconnected 
flux, and region II represents reconnected flux. The two regions I are sepa- 
rated by the very narrow reconnection current layer. Plasma from regions I 
enters the reconnection layer and gets accelerated along the layer, finally en- 
tering the separatrix region, also between regions I and II. In general, both 
the reconnection layer and the separatrix region require resistive treatment. 

The plasma entering the separatrix from the reconnection layer is trav- 
eling at nearly the Alfven speed. It crashes into the plasma at rest on the 
separatrix lines, that has not passed through the reconnection layer, but 
which got there by direct E x B motion across the lines as the position of 
the separatrix changed by reconnection. This crash generates considerable 
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Separatrix : 
region 

Figure 1: The global geometry. 

heat, and hence back pressure on the reconnection layer. However, the sep- 
aratrix region is continually in transition since different plasma occupies it 
as the reconnection proceeds so that the heated plasma is moved into the 
reconnected region. Thus, the plasma encountered by the outflowing recon- 
nected plasma is continually refreshed and can always be taken initially at 
rest. Therefore, there is a time delay before the back reaction sets in. In 
our paper, we attempt to model this dynamical behavior of the separatrix 
plasma as accurately as possible. 

In the limit S —>■ oo the reconnection rate is slow compared with the 
Alfven time ta = L/Va, which allows one to break the whole problem into 
the global problem and the local problem. In a previous paper Ref. 0, we 
argued that on the global scale (i.e., on the scale of order the half-length of 
the layer L) the time evolution of the reconnecting system can be described 
as a sequence of magnetostatic equilibria. In paper || we explained that the 
role of the global solution is to give the general geometry of the reconnecting 
system, the position and the length of the reconnection layer and of the 
separatrix, and the boundary conditions for the local problem (which, in turn, 
determines the reconnection rate). These boundary conditions are expressed 
in terms of the outside magnetic field B y fi(y), where y is the direction along 
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the layer. In particular, B y fi(y) provides the characteristic global scales: the 
half-length of the layer L, defined as the point where By^iy) has minimum, 
and the global Alfven speed, defined as Va = B y>0 (0) / 'y/ 'An p. It is important 
to understand that the global solution is essentially independent of the local 
reconnection physics. 

In this paper we study the local problem concerning the reconnection 
layer itself. Our main goal here is to determine the internal structure of a 
steady state reconnection current layer (i.e., to find the 2D profiles of plasma 
velocity and magnetic field), and the reconnection rate represented by the 
(uniform) electric field E. Our most important result is that, in the case of 
uniform resistivity, there is a unique stable steady-state solution, and this 
is essentially the Sweet-Parker solution. We show that Petschek's solution 
(which is not unique and supposedly encompasses a variety of reconnection 
solutions from that of Sweet and Parker to a solution reconnecting at almost 
the Alfven speed) actually relaxes to the single Sweet-Parker solution. 

First, in Section |2| we discuss the assumptions of our physical model of the 
layer in some detail. Then, in Section | we introduce the rescaled equations 
representing the mathematical model of our problem. In Section |] we present 
our numerical simulations. And, finally, in Section |5| we give our conclusions. 

2 Physical Model 

By the local problem we mean the analysis of the internal structure of the 
reconnection layer and the separatrix layer, which is necessary for the de- 
termination of the reconnection rate. Since in ideal MHD these layers are 
current sheets of zero thickness, resolving their inner structure cannot be 
done in this ideal framework and requires the addition of some new nonideal 
physics. 

Historically, most important, and conceptually most interesting for the 
purpose of resolving the current layer, is the inclusion of small resistivity. 
The model in which the only nonideal effect is that of the resistivity appears 
to be the simplest model with the minimal required complexity needed to 
resolve the current density singularity of an ideal MHD solution. 

Thus, in this analysis we assume that the only new nonideal physical 
process is small constant and uniform resistivity r\ (and perhaps viscosity u; 
see discussion below). 

The inclusion of the small resistivity means the introduction of a new di- 
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mensionless small parameter associated with it, namely the inverse Lundquist 
number S' 1 = tj/LVa <C 1, which is considered the primary small parameter 
of the problem. This means that we are interested in studying the case of 
S ^ oo, and we want to determine how the parameters of the layer, such as 
its thickness and the reconnection rate, scale with S in the leading order. As 
we shall see in the next section, making use of this small parameter helps to 
simplify the problem significantly while keeping all essential features, such 
as the two-dimensional nature of the problem, intact. 

The next assumption we make concerns a steady state. The steady-state 
condition means that parameters of the current layer as well as the boundary 
conditions change very slowly compared with the global Alfven time, which 
is the characteristic time spent by a fluid element inside the reconnection 
layer. A very important consequence of the steady-state condition in 2D 
geometry is that, due to the Maxwell equation V x E = — (1/c) dH/dt, the 
z component of the electric field is uniform: E z (x, y) = E = const. 

As for the plasma viscosity, it does not seem to be necessary to include 
it, because viscosity, unlike resistivity, does not play any role in the actual 
breaking the lines of force. However, we always include a small constant and 
uniform viscosity v for two reasons. First, in our numerical simulations we 
include it for numerical stability. The second and more important reason 
is that the consideration of very small (even compared with the resistiv- 
ity) viscosity is useful to correctly understand some important features of 
the magnetic field configuration at the very center of the current layer (see 
Ref. I). 

Next, most of the classical models of reconnection, including both Sweet- 
Parker and Petschek, assume that the outside merging magnetic field is uni- 
form. But this assumption actually prohibits one from formulating the prob- 
lem in a mathematically complete and consistent manner, because the down- 
stream boundary conditions for the flow cannot be correctly specified. In this 
quasi-one-dimensional framework, there is no natural end of the layer; in par- 
ticular, there is no way to define the global scale L. This, in turn, makes all 
attempts to get some scalings for the reconnection rate with the Lundquist 
number essentially meaningless, since the definition of the Lundquist number 
involves L. 

Now, in our paper a generic and more or less representative variation of 
the outside magnetic field B yt o(y) along the layer and along the separatrix 
[where it is called B St i(l)] is included as an integral part of the problem. In 
particular, the global scale (the half-width L of the layer) is defined naturally 
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as the distance along the midplane from the center of the current layer to 
one of the two endpoints — i.e., the points where the outside magnetic field 
goes through a minimum and where the separatrices branch off the midplane 
x = 0. In fact, this global scale is the characteristic scale for the function 
B y fi{y). Since this function is determined by the global ideal MHD solution 
(see Ref. ||), the scale L is, by definition, independent of the physics of the 
resistive layer. 

Thus, the nonuniformity of the outside magnetic field along the layer 
makes the problem essentially two-dimensional (rather than one-dimensional). 
One practically important consequence of this fact is that the problem be- 
comes much more complicated mathematically, so that one has to abandon 
any hope for a nice analytical solution and to resort to numerical simulation 
instead. 

Thus, the physical model of the reconnection layer that we are going to 
use in this paper for treating the local problem, can be summarized as 

two-dimensional, steady-state, incompressible, resistive MHD with 
constant and uniform resistivity (and perhaps viscosity) in the 
limit of very large Lundquist number. 

Perfect mirror symmetry is assumed with respect to both the x axis and 
the y axis.[] 

We call this model the canonical reconnection layer model. 

When considering the local problem in this model, we use the global scale 
in the y direction (along the layer), and the local scale in the x direction 
(across the layer). The outside magnetic field B y> o(y) determined by the 
global solution here plays the role of a boundary condition at x — > oo. 

Finally, although in this section we talked only about the reconnection 
layer itself, the same physical model applies also to the separatrix layer. 

3 System of Rescaled Equations 

In accordance with our physical mode, we can now write down the set of two- 
dimensional steady-state fluid equations for our system. These equations are: 

1 Duc to this symmetry, v y and B x are even and v x and B y are odd with respect to 
x = 0, and v x and B y are even and v y and B x are odd with respect to y = 0. Thus, v x = 
along x = and v y = along y = 0, which means that the two axes of symmetry are 
stream lines. 
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(i) The incompressibility condition: 



V .v = ^ + ^ = 0. (1) 

ox oy 

(ii) The z component of Ohm's law: 

i]j z = E + v x B y - v y B x , (2) 

where E = E z = const. 

(hi) The equation of motion (with the viscosity): 

v ■ Vv = - Vp + [j z x B] + zA7 2 v, (3) 

where the density is set to one. 

Now we take the crucial step in our analysis. We note that the reconnec- 
tion problem is fundamentally a boundary layer problem, with S~ l being the 
small parameter. This allows us to simplify our MHD equations by perform- 
ing a rescaling procedure^L^ inside the reconnection layer, to make rescaled 



resistivity equal to unity. This can be done in a natural way if one rescales 
the distances and the fields in the y direction to the corresponding global val- 
ues (i.e., the length of the layer L, the outside magnetic field just above the 
center of the layer B y ^(0), and the corresponding Alfven speed Va), while 
rescaling the distances in the x direction and the x components of the velocity 
and magnetic field to the corresponding local values: 

y x 

x, 



L " 5o 



V A ' Vyi VaSo/L 
B,,i B 



B V M B y>0 (0)8o/L 

Pi „ TT^Tr f it E > 



Bl (0)/4n " B yfi (0)V A 5 /L 

where = LS' 1 ^ 2 is the Sweet-Parker thickness of the current layer. Thus, 
one can see that the small scale Sq emerges naturally as the thickness of the 
resistive boundary layer. 
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The viscosity v is now rescaled as u/rj — > z/. We assume that it is at 
least as small as the resistivity [which means v = 0(rj) or less], and most of 
the time (see Ref. ||) we will be interested in the case of vanishing viscosity 
v — > (which means v <C rf). 

Now, all of the rescaled dimensionless quantities (v x , v y , B x , B y ) are 
generally of order 1. Using the small parameter 5q/L = S^ 1 ^ 2 <C 1, one can 
simplify the equations by writing them down in the leading nontrivial order 
in Sq/L. This way one neglects all unimportant corrections, keeping only the 
essential terms. 

First, the incompressibility condition is written in rescaled quantities in 
exactly the same way as in the unrescaled quantities: 

^ + ^ = 0. (5) 
ox oy 

The z component of the steady-state Ohm's law can be written as 

dB 

E = - v x B y + v y B x , (6) 

where the first term on the RHS is the resistive term. 

Next, consider the equation of motion. Since all the velocities in the 
x direction are small compared with the Alfven speed, the inertial terms in 
the x component of the equation of motion (3) are small, and this equation 
just gives one the pressure balance across the current sheet: 

= 0, (7) 

which allows one to determine the pressure p in terms of B y (x,y) once the 
pressure and the magnetic field B y< o(y) outside the reconnection layer are 
known. It is customary to set the pressure outside the layer to zero, so that 

Bl (y) Bfa,y) 

p(x,y) = ^ ^ — • (8) 

(Of course, one can add a large constant to the pressure to enforce incom- 
pressibility. Since this constant cancels out of all our equations, we ignore 
it.) 
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Finally, one has the y component of the equation of motion, with acceler- 
ation provided both by the pressure gradient force, and by magnetic forces, 
and with the viscous force: 

„ dp _ dB v d 2 v v 

ay ox ox 1 

It is interesting to note that almost everywhere in the layer the magnetic 
force j z B x = B x d x B y is actually just the y component of the magnetic pres- 
sure gradient force V±B 2 /8n = VB 2 /8tt — bb • VB 2 /8tt where b is the unit 
vector along B. The magnetic tension force B 2 h ■ Vb/47r is higher order 
everywhere except very close to the midplane (in a thin region where the 
unrescaled fields satisfy B y < B x ). 

We believe that this rescaling procedure captures all the important dy- 
namical features of the reconnection process. 



4 Numerical Simulations 

In order to find the steady-state solution for the system, we designed a time 
dependent resistive MHD code for the main reconnection layer. This main 
code was supplemented by another code describing the separatrix, which is 
needed to provide the nearly correct downstream boundary conditions for 
the main code. 

In Section [4.1| we describe the main code for the reconnection layer to- 



gether with the boundary conditions. In Section |4.2| we discuss the different 
choices for the initial conditions used in the simulation. In Section |4l| we 
present the model for the separatrix region that we have used in order to get 
the downstream boundary conditions for the main layer. In Section [O we 



report the results of our numerical simulations. More details can be found 



in Ref. 1C . 



4.1 Numerical Scheme and Boundary Conditions 

In order to approach the steady-state solution described by the system of 
rescaled equations (5), (6), (8), and (9), we followed the true time evolution 
of the system, starting with some initial conditions that will be described in 



Section 4.2 
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The time evolution is described by the following two dynamic equations 
for the two dynamic variables ^ and v y : 

^ = _ v . (v *) + _ + ^_j, (10) 

and 



+ V .(B^) + ^ (11) 



These equations are written (using V-v = 0) in conservative form (i.e., in the 
form of conservation laws), which is preferable for numerical computations. 
Small artificial resistivity r\ y and viscosity v y , acting in the y direction, are 
added to provide numerical stability. Because they are small, these terms do 
not change the solution noticeably, as was verified in the runs. The natural 
unit of time in our simulations is the global Alfven time T4 = L/Va- 

Once v y and \l/ are known everywhere at a new time step, one can find 
all other variables. Namely, B x and B y are given by the derivatives of \I/: 



B * = -dj> By = a± ' (12) 



and v x is obtained from the incompressibility condition by integrating dv y / dy 
in the x direction starting from the midplane, where v x = because of the 
symmetry: 

v x (x,y) = - -—^dx. (13) 
Jo ay 

We used the finite-difference method with centered derivatives (providing 
second order accuracy) in both the x and the y directions. 

The time derivatives were one-sided. Our scheme was explicit in the 
y direction, but in the x direction the resistive term in Ohm's law (d 2 ^/dx 2 ) 
was treated implicitly, while all other terms were treated explicitly. This 
enabled us to speed up the computations. 

We conducted the simulations on a rectangular grid (i m x j m ). Because 
of the symmetry, we considered only one quadrant (see Fig. |^). 

Now let us discuss the boundary conditions. There are four boundaries 
in the system: the upstream (or upper) boundary x = X\ im , < y < yu m , 
the downstream (or right) boundary y = y\ im , < x < X\ im , and the two 
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Figure 2: Computational box used in the numerical simulation. 



boundaries formed by the axes of symmetry: the lower boundary x — 0, 
< y < yii m , and the left boundary y = 0, < x < x hm . The flow 
enters through the upper boundary and leaves through the right boundary, 
so that there is no flow of plasma through the left and the lower boundaries. 
While the boundary conditions at the left and the lower boundaries come 
from simple symmetry conditions, the boundary conditions on the upper 
boundary and especially on the right boundary are more complicated, as we 
shall discuss below. 

On the upper boundary, the boundary conditions come from matching 
with the ideal MHD solution in region I above the reconnection layer. In 
rescaled quantities, this matching should be done at x — > oo. But since in a 
numerical simulation it is not possible to place a boundary of the computation 
box at infinity, we place it at some sufficiently large xu m 1. The typical 
values of x\ im in our simulations were rri im ~ 4 — 8. 

From the ideal solution in region I we know that, as x — > oo, v y — > 
(meaning v y <C Va) and B y — > B y (y), which is prescribed. It turns out that, 
since the upper boundary is placed not at infinity but at some finite (although 
large) xu m , it is better, for numerical reasons, to choose dv y /dx(xu m ) = 
instead of v y (xu m ) = 0. 
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As for the magnetic field, we just set d^$! / dx(xu m ) = B y0 (y). The function 
B y fi(y) depends on the external boundary and the thermodynamics of the 
global solution as well as the amount of reconnected flux. Rather than make 
use of a specific global solution we choose a generic one since the physics in 
the reconnection layer should be the same. However, near y = L (the end of 
the reconnection layer which terminates in a cusp geometry we must employ 
the square- root behavior derived in Ref. ||, which any global solution must 
have. Thus, in our numerical simulations we typically took 

B v ,o(y) = B + (1 - B )yjl-y*, (14) 

consistent with the cusp solution (see Ref. ||). The value B of the outside 
magnetic field B y $ at the endpoint y — 1 typically was taken to be Bq = 0.2 
or B = 0.3 [the magnetic field is normalized so that B yj0 (0) = 1]. 

This choice of boundary conditions worked well in our simulations. In 
particular, the behavior of the solution near the upper boundary was smooth, 
and the solution deep inside the reconnection layer did not depend on the 
exact position of the upper boundary. 

At the lower boundary x = 0, the boundary conditions come naturally 
from the requirement that both v y and \l/ be symmetric with respect to the 
midplane. Thus, 

^ = and B y = ^- = (x = 0), (15) 

and also v x (0,y) = 0. 

Here, however, we would like to make one remark. The boundary con- 
dition for dv y /dx at x = is needed only when one includes the viscous 
term vd 2 v y jdx 2 in the equation of motion. If one does not keep this term, 
then this equation contains only first derivatives of v y in the x direction, so 
one needs only one condition, which can be set at the upper boundary. In 
our simulations, however, we always include viscosity (usually small, but not 
zero), both for numerical reasons and in order to resolve the behavior near 
the midplane (see Ref. ||). 

The boundary conditions at the left boundary are similar to those on the 
lower boundary. They follow from symmetry and are rather straightforward: 

v y (x,0) = 0, and B x (x,0) = 0. (16) 

Now consider the right (or downstream) boundary. We chose it at some 
point y = y\ im close to the endpoint (typically y\ im = 0.9 — 1.0). Then, one 
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needs to specify the downstream boundary conditions on this boundary. This 
boundary is, in fact, the interface between the main layer and the endpoint 
region and the separatrix. The boundary conditions should describe the 
effect of the separatrix reacting back on the main layer, in particular the back 
pressure. We are not aware of any previous numerical or theoretical studies 
in which the role of the back pressure has been adequately investigated. The 
problem of how to set the boundary conditions on this boundary is rather 
nontrivial and its discussion is postponed until Section fO|. 

To summarize, the advantages of our approach to numerical simulation 
of the reconnection layer are the following: 

1) First, the use of rescaled equations takes us directly into the realm of 
S — > oo. 

2) Second, this is an essentially 2D (rather than ID) code that uses a 
realistic variation of the outside magnetic field along the layer. The position 
of the endpoint is clearly defined in terms of the function B y)0 (y). We do not 
assume B y fl(y) = const as many people do. 

3) We obtain the steady-state solution by following the true time evolu- 
tion, and the rescaled equations are such that we do not give the boundary 
conditions for the incoming flow velocity at the upstream boundary. This is 
because the physical v x is small and its evolution is not determined from a dy- 
namic equation of motion. Instead, the equation of motion in the x direction 
simply degenerates into the vertical pressure balance, and v x is determined 
from the incompressibility condition. Thus, we do not specify v x (xn m , y) as a 
boundary condition, which means that we do not prescribe the reconnection 
ratel The system itself determines what the reconnection rate should be! 
This is really a very important point. 

The fact that we rescaled the x coordinate using the Sweet-Parker scaling 
does not actually mean that we prescribe the Sweet-Parker reconnection rate. 
If the system wants to go at a faster rate, then it would try to develop some 
new characteristic structures, extending beyond Xn m , and we should be able 
to see it. 

4.2 Initial Conditions 

We have performed several runs with different initial conditions consistent 
with our boundary conditions. 

In some cases we started with a configuration qualitatively resembling the 
Sweet-Parker reconnection layer (see Figs. |3] and f|). These initial conditions 
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can be written in the following analytical form: 



B y (x,y) = B yfi (y) tanh(x), (17) 
V(x,y) = V (0,y) + B y , (y) log cosh x, (18) 

B x (x,y) = - — , (19) 
dy 

where we took the outside magnetic field B y (y) in the form (14), and the 
variation of the magnetic flux on the midplane as 

*(0,*/) = -^/ 2 - (20) 
The initial velocity was taken in the form 

v x (x, y) = - E ° tanh(x), (21) 
By,o{y) 

and 

v y (x,y) = -J q ~^dy, (22) 

where the initial outside electric field Eq varied but typically was of order 
one. We call the initial conditions described by the set of Eqs. (4.17)-(4.22) 
the Sweet-Parker-like initial conditions. 

In other runs we wanted to see whether the system would want to go with 
the faster reconnection rate and whether it would develop the Petschek-like 
structures. 

(It should be observed that Petschek, in his original paper, did not spec- 
ify a unique reconnection rate. Instead his reconnection rate depended on 
the length of the dissipation region y*, and his reconnection velocity was 
^rec = (\j L/y^VA/yfS. For example, the choice y* = L gives the Sweet- 
Parker rate. The usually quoted Petschek reconnection velocity, fp e tschek ~ 
Vy^/logS, was obtained by him when he took the minimum possible value 
for that would not seriously perturb the global solution. Furthermore, it 
is possible to include Petschek's shocks inside the boundary layer as indeed 
he did.) 

We conjectured that if the system wants to go at a faster reconnection 
rate than that of Sweet-Parker, then it would at least be able to go at a rate 
twice as fast(ie., y* = L/4). Therefore, we carried out several runs where at 
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Figure 3: Contour plot of the magnetic field for Sweet-Parker-like initial 
conditions. 
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Figure 4: Plots of j(x) for different values of y at t — for Sweet-Parker-like 
initial conditions. 
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t = we set up a Petschek-like structure (see Figs. ||, and [7]) described by 
the following expressions: 



*(0,y) = -— log cosh (2; ), 



B y (x,y) = B y!0 (y) 



tanh(zi) + tanh^) 



(23) 
(24) 



V(x,y) = ^(0,y) + %^ pogcosh(z 1 )+logcosh(z 2 )-21ogcosh(zb)], (25) 



2M 



v y (x,y) 



tanh(£i) — tanh(z 2 ) 



v x (x,y) = —M 



tanh(zi) + tanh(z 2 ) 



(26) 
(27) 
(28) 



where Zq = M 2 y, Z\ = Mx + M 2 y, and Z2 = Mx — M 2 y, and where the pa- 
rameter M corresponds to the initial reconnection rate in terms of the Sweet- 
Parker reconnection rate. It describes how well pronounced the Petschek-like 



structure is (M would correspond to yL/y* in Petschek's notation). Typi- 
cally M was chosen to be 2 or 3. We call the initial conditions described by 
the set of Eqs. (23)-(28) the Petschek-like initial conditions. 



4.3 The Downstream Boundary Conditions and the 
Model for the Separatrix Region 

The downstream boundary y = yu m is the interface between the main layer 
and the separatrix region. The boundary conditions at this boundary cannot 
be given in a simple closed form.f] Instead, they require matching with the so- 
lution in the separatrix region, which itself is just as complicated as the main 
layer. Therefore, we have developed a supplemental numerical procedure for 
the separatrix region. 

2 For example, we have tried to use the so-called free-flow boundary conditions: 
d 2 ^/dy 2 = 0, d 2 v y /dy 2 = 0. The steady-state solution exists and is reached within 
several Alfven times. However, the solution in the bulk of the main layer strongly depends 
on the position of both the upstream and downstream boundaries, which is physically not 
acceptable. 
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Figure 5: Contour plot of magnetic field at t — for the Petschek-like initial 
conditions. 




Figure 6: The current density j(x,y) at t — for Petschek-like initial condi- 
tions (all four quadrants are shown for clarity). 
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Figure 7: Plots of j(x) for different values of y at i = for the Petschek-like 
initial conditions. 

To make the situation more tractable, we have adopted simplified physical 
model for the separatrix — namely, a model in which the resistive and the 
viscous terms are omitted, and hence the magnetic field is frozen into the 
plasma. Even though this model does not describe the separatrix completely 
accurately, it should give one a qualitatively correct picture of the influence of 
the separatrix region back on the main layer, and thus sufficiently reasonable 
downstream boundary conditions for the main layer. In particular, our model 
includes the effects of the back pressure that the separatrix exerts on the main 
layer. 

We follow one given field line as it moves across the separatrix of length L s 
as described in the introduction. This field line is described in terms of two 
functions: the magnetic field B(l, t), and the parallel velocity v(l, t) = v»(l, t). 
Here I is the length measured along the field line, starting from the boundary 
between the separatrix and the main layer (I = or y = yn m ), and ending at 
the reflection point A at the top of Fig. [I] (I = L s ). 

Now, the time evolution in this one-dimensional problem corresponds 
to the perpendicular motion of the field line through the separatrix region. 
The new time variable t represents the relative position of the field line in 
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real 2D space with respect to other field lines, and there is a one-to-one 
correspondence between the time t and the x position of the footpoint (the 
point where the field line intersects the boundary y = yu m ) of this field 
line. This correspondence between t and x is given by the x component of 
the velocity of the field line, v xs = —E/By(x,yu m ) < (x decreases as t 
increases), and by the initial condition that x = Xu m at t — 0. We stop 
the simulation of the separatrix at the moment t when x(t) becomes zero - 
i.e., when the midplane is reached. A substantial advantage of this model is 
that one actually does not need to know where the field line is in real 2D 
space! The only thing one does need to know about the position of the field 
line in order to set up the downstream boundary conditions for the main 
layer is where it connects to the main layer - - i.e., what is the value of 
the coordinate x at the point where the field line intersects the downstream 
boundary y = y\ [m (i.e., the footpoint). 

The boundary conditions on the upper boundary of the separatrix are 
v (I) = and B{1) = B sl (t), where B s i(t) is the field outside the separatrix 
in the upstream region given by the global ideal MHD solution, as explained 
in Ref. Hi- In our simulation we took this outside field in the form 



B,Al) 



\ 











I 1 - 2L,j 



(29) 



Now these conditions become the initial conditions at t = for our ID 
(plus time) problem concerning this one field line. 

As in the main layer, we use here an explicit (in /) code with centered 
differences for convective derivatives in the / direction. The boundary con- 
ditions at the right boundary I = L s (the reflection point A on Fig. [I]) come 
from the condition of symmetry with respect to this point (I = L s ): 

v(L SJ t) = } _(L s ,t)=0. (30) 

The other pair of boundary conditions is given at the left (incoming) 
boundary I = and is provided by the main layer in terms of the values of 
v y (x,y\i m ) and B v (x,yu m ). Assuming that the steady state in the main layer 
has already been achieved, these boundary conditions are constant in time in 
the laboratory frame, but in the frame moving along together with the field 
line, they are now time-dependent. 
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Since at the boundary y = y lim the layer is still essentially straight, and 
the field lines (except those few very close to the midplane, which are ignored) 
are also almost straight, then the absolute value of the magnetic field B is 
almost exactly equal to B y and the parallel velocity v is almost exactly equal 
to v y . Therefore, the matching conditions at this boundary are B y = B, 



v y 



As for the downstream (or lower) boundary, one really does not need to 
set any conditions there because the x coordinate in our model corresponds 
to time and the new equations have only first-order time derivatives. 

Now let us derive the differential equations for this model of the separa- 
trix. The main idea in this derivation is that the convective term vi • Vj_ is 
replaced by the time derivative d/dt. The natural unit of time in this model 
is the Alfven time, the natural unit of distance in the I direction is the global 
scale L, and that of the parallel velocity is the Alfven velocity Va- 

First of all, one can apply this model only to the part of the separatrix 
where the field lines are not very strongly curved (i.e., where the radius of 
curvature of the field lines is of order L) , the magnetic field itself is sufficiently 
strong (of order B Q or at least B 3> B 1 = E/Va), and the perpendicular 
velocity v± (which in ideal MHD is equal to E/B) is small compared with 
the parallel velocity v\\ = 0(Va)- Thus, the model is bound to fail in the 
very small region near the midplane where the real physical B y becomes 
comparable with B x ^\ 

In the system of reference moving together with the field line in the direc- 
tion perpendicular to B, the (parallel component of) the magnetic induction 
equation becomes 

B=^ = B±v-4b, (31) 



dt dl dl 
and the parallel equation of motion becomes 



dv dv d 

at = ~ v ~m ~ dl 



B 2 sl (l) B*(l,t) 



(32) 



Eq. (32) is equivalent to Eqs. (7) and (11). Physically, the force parallel to 
B is the pressure gradient force which can be expressed in terms of B by 



3 This presumably occurs at some infinitesimally short distance 5i from the midplane, 
where 5\/8q ~ Sq/L <C 1. Since one considers the limit S — > oo, and hence So/L — > 0, 
this region shrinks to zero, and one should not be concerned too much about it in our 
simulation. 



20 



Eq. (7). It is interesting to note that, while the fluid in the initial problem 
was incompressible, in this ID problem the motion is effectively not incom- 
pressible (dv/dl 7^ 0), and the thickness of the flux tube (~ l/B) plays the 
role of density. 

Now let us see how the back pressure from the separatrix acts on the 
main layer. 

As the field line moves, the incoming velocity v(l = 0,t) increases from 
almost zero at t — to about the Alfven velocity, and the incoming magnetic 
field B{1 = 0,t) drops from Bq to zero at the midplane. Therefore, there 
is a point x = x c somewhere in the middle of the left boundary of the 
separatrix where the incoming flow becomes locally super- Alfvenic (v > B). 
The propagation of information with respect to the fluid in our model occurs 
at the local Alfven speed by the means of two characteristics dl/dt — v ± B. 
For x > x c (v < B) , one characteristic goes from the left boundary I — (y — 
2/iim) towards the right boundary, while the other characteristic goes from the 
upper boundary x = X\ im towards the left boundary, carrying information 
about the pressure and the flow in the separatrix to the main layer. This 
means that, for x > x c , the layer "feels" the effect of the separatrix in the 
form of the back pressure coming from the previously undisturbed fluid in 
the separatrix region. This back pressure is found to have some stabilizing 
effect, making the solution deep inside the main layer independent of x hm and 
almost independent of the position y iim of the downstream boundary (where 
the separatrix and the main layer solutions are matched). 

After the field line crosses the point x = x c , the flow becomes super- 
Alfvenic and the characteristic dl/dt = v — B > is deviated away from the 
left boundary back into the separatrix region. Now both characteristics come 
out of the left boundary and there is no propagation of information from the 
separatrix region into the main layer. This also means that when the field 
line reaches the vicinity of the midplane x = and the description of the 
separatrix region as a region of almost straight field lines fails, it does not 
matter much, because the flow at this point is very strongly super-Alfvenic 
and the main layer does not feel what happens downstream. We would like 
to remark that in our simulations we have not observed any shock formation 
in the separatrix. 

Of course, this picture of the flow of information is valid only in this ideal 
MHD model. In the real situation with resistivity, there is always propagation 
of information upstream due to the resistive diffusion. 

To summarize, the model of the separatrix region presented here has 
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two main drawbacks. First, it describes an ideal MHD separatrix. Sec- 
ond, it assumes that the boundary conditions on the incoming boundary 
I — (y — 2/ii m ) are stationary (i.e., that the main layer has already reached 
the steady state). (This second objection can be removed if the main layer 
is evolved with right hand boundary conditions changing in step with the 
separatrix solution and vice versa, but this would involve a longer simula- 
tion.) In addition, the model is valid only when the field lines are essentially 
straight, as discussed above (the region where this assumption breaks down 
is discussed in Ref. |10|). Despite all that, however, we feel that this model 



provides a qualitatively correct picture of the dynamics in the separatrix re- 
gion in the steady state, and thus gives us sufficiently reasonable downstream 
boundary conditions for the main layer. 

4.4 Results of the Simulations 

Now let us present and discuss the results of our numerical simulations. 

We found that, after a transient period of a few Alfven times, the system 
reaches a steady state that is independent of the initial configuration. 

In particular, when we start with a Petschek-like initial configuration 
(described in Section |4~2j ), the high velocity flow rapidly sweeps away the 
transverse magnetic field B x (see Fig. []). This is important, because, for a 
Petschek-like configuration to exist, the transverse component of the mag- 
netic field on the midplane, B x (0,y), must be large enough to be able to 
sustain the Petschek shocks in the field reversal region. For this to hap- 
pen, has to rise rapidly with y inside a very short diffusion region, B x (0,y) 
y < y* <C L (in the case M = 2, presented in Fig. [| y* — L/4), to reach a 
certain large value (B x = 2 for M = 2) for y* <C y < L. While the transverse 
magnetic flux is being swept away by the plasma flow, it is being regener- 
ated by the merging of the B y field, but only at a certain rate and only on 
a global scale in the y-direction, related to the nonuniformity of the outside 
magnetic field B y ^(y), as discussed in Refs. and |T0| . As a result, the ini- 



tial Petschek-like structure is destroyed, and the inflow of the magnetic flux 
through the upper boundary drops in a fraction of one Alfven time. Then, 
after a transient period, the system reaches a steady state consistent with 
the Sweet-Parker model. 

In general, the transient period typically lasts a few Alfven times, during 
which the incoming electric field can oscillate around its final steady-state 
value. These oscillations have a period of order ta, and a decay time also 
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Figure 8: The time evolution of the variation of the transverse magnetic field 
B x (0,y) along the midplane x = for the Petschek-like initial conditions. 

of the same order. After several ta, the electric field becomes constant and 
uniform throughout the computational domain and the system approaches 
steady state. We terminate our simulations typically after 5 or 10 Alfven 
times. The time evolution of the incoming electric field [i.e., E(x = xu m , 0)(t)] 
for different choices of the initial conditions is represented in Figs. |9| and [K| 
We believe that the fact that we rescaled x using the Sweet-Parker scaling 
does not mean that we prescribe the Sweet-Parker reconnection rate. Indeed, 
if the reconnecting system wanted to evolve towards Petschek's fast recon- 
nection, it would then try to develop some new characteristic structures, e.g., 
Petschek-like shocks, which we would be able to see. Note that, if Petschek 
is correct, then there should be a range of reconnection rates including those 
equal to any finite factor greater than one times the Sweet-Parker rate E$p. 
However, in our simulations we have demonstrated that there is only one 
stable solution and that it corresponds to E = E$p. In this sense we have 
demonstrated that Petschek must be wrong since reconnection can not even 
go a factor of two faster than Sweet-Parker, let alone almost the entire factor 
oiy/S. There seems no alternative to the conclusion that fast reconnection 
is impossible. 
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Figure 9: The time evolution of the electric field E at the point x = xu 
y = yiim for the Sweet-Parker-like initial conditions with E = 1.0. 



E(t) 




Figure 10: The time evolution of the electric field E at the point x = xu m , 
y = y\im for the Petschek-like initial conditions with M = 2 (corresponding 
to E = 2.0). 



24 



It is interesting that in Petschek's original paper the length of the 
central diffusion region y* is an undetermined parameter, and the reconnec- 
tion velocity t> rec depends on this parameter as VA{L/y*) 2 / \/~S . If is taken 
as small as possible then Petschek finds that v rec ~ Va/ log(S). However, 

should be determined instead by balancing the generation of the trans- 
verse field B x against its loss by the Alfvenic flow (it should be remarked 
that Petschek did not discuss the origin of this transverse field in his pa- 
per). As we discussed above, this balance yields ~ L, with the resulting 
unique rate equal to that of Sweet-Parker. This results are borne out by our 
time-dependent numerical simulations. 

The final steady-state configuration represents the Sweet-Parker recon- 
nection layer. This means that all the plasma parameters are of order one in 
the rescaled coordinates, and change on a scale of order <5 in the x direction 
and on a global scale L in the y direction. For our choice of boundary con- 
ditions [i.e., of the outside magnetic field B y<0 (y)], the reconnection rate in 
the steady state was E = 1.0-Bgp, where E SP = r/ 1 / 2 V^-B^o (0) is the typical 
Sweet-Parker reconnection rate. 

The use of the time-dependent equations allows us not only to find the 
steady-state solution, but also to draw some conclusions about its stability. 
The Sweet-Parker solution was found to be stable and robust: it did not 
depend on the positions of the boundaries xn m , yu m or on the small artificial 
resistivity and viscosity. Moreover, we found that it is fairly insensitive even 
with respect to the choice of the parameters describing the outside magnetic 
field, such as B = B y> o(l) (thus, E varies by about 10% as B changes from 
to 0.3). 

The steady-state Sweet-Parker solution is represented in Figs. O, O, [13], 



1% and [15|. This solution corresponds to the following set of the parameters 
describing the boundary conditions: x\ im = 5.0, yn m = 1.0, B = 0.3, L s = 
2.0, C = 0.957 [see Eqs. (14) and (29)]. In this particular run the values of 
the perpendicular viscosity u, and of the (artificial) resistivity and viscosity 
acting in the y direction were v = 0.02, rj y = v y = 0.01. 

There are several things that should be noted about this solution: 
As can be seen from Fig. 13, the current density as a function of x at 
any given value of y peaks on the midplane x — 0. This makes the solution 
qualitatively different from a Petschek-like configuration (see Fig. [j]) in which 
the current density is concentrated in a shock-like structure off the midplane. 

At any given y, the current density j(x) rapidly goes to zero as x — > oo, 
and B y (x, y) — > B y (y) monotonically, which means that there is no flux pile 
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Figure 11: Contour plot of the steady-state magnetic field in the reconnection 
layer. 



up in front of the layer. The velocity v y does not quite go to zero at the 
upper boundary, but its value at x — xn m is small and goes to zero as the 
artificial r] y , u y — > 0. 

Next, the solution in the layer shows an essentially linear rise of v y and 
B x along the midplane x = (see Figs. [16| and [IT]). The linear behavior of 
B x {x = 0,y) near y = 0, contrary to the cubic behavior predicted by Priest 
and Cowley |]TT |, was explained in Ref. |§ (together with the nonanalytic 
behavior of the solution near the midplane in the limit v — > 0). 

Finally, B x exhibit a sharp change near the downstream boundary y = 
y\ im , as can be seen in Figs. [L7[ This change is due to the fact that in the 
separatrix region we neglect the resistive term ^ xx , which is in fact finite. 
That is, the (perpendicular) resistivity effectively has a discontinuity across 
V = Vum- v(y < yiim) = l and 7](y > y hm ) = 0. This discontinuity in the 
equations also shows up in the solution, but it is smoothed out over some 
vicinity of yn m by the artificial resistivity and viscosity in the y direction. 
As these r] y , v y go to zero, the region of the rapid change near y\ im becomes 
smaller and smaller. 

Let us add a few remarks about the role of the separatrix region in our 
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Figure 13: Plots of j as a function of x for several different values of y in the 
steady state. 
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Figure 14: Plots of B y as a function of x for several different values of y in 
the steady state. 
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Figure 15: Plots of V y as a function of x for several different values of y in 
the steady state. 
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Figure 16: Plots of v y as a function of y along the midplane x = in the 
steady state. 




Figure 17: Plots of B x as a function of y along the midplane x 
steady state. 
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simulations. 

First, it appears that the destruction of the initially-set- up Petschek-like 
configuration and its conversion to the Sweet-Parker-like layer is so robust 
and happens so fast, that it is determined by the dynamics in the main 
layer and by its interaction with the upstream boundary conditions [i.e., 
with the scale of nonuniformity of the outside magnetic field B yt0 (y)]. As a 
consequence, it has nothing to do with the downstream boundary conditions 
(i.e., with the separatrix region). Therefore, the fact that our model of the 
separatrix does not describe the separatrix completely accurately seems to be 
unimportant, as far as the instability of the Petschek solution is concerned. 

Thus, we believe that the separatrix region, while providing physically 
reasonable downstream boundary conditions for the main layer problem, does 
not really have a strong effect on the principal result that the stable steady- 
state solution is the Sweet-Parker layer with the Sweet-Parker reconnection 
rate. Still, we have to point out that for the solution of the problem to be 
really complete, one needs to build a better model of the separatrix region. 
Such a model would include real time dependence and resistivity, and also 
would treat the very near vicinity of the endpoint of the reconnection layer. 
A proper consideration of the endpoint region cannot be done in rescaled 
variables, and a further rescaling of variables and matching is needed. 

5 Conclusions 

To summarize, in this paper we present a definite systematic solution to a par- 
ticular clear-cut, mathematically consistent problem concerning the internal 
structure of the reconnection layer within the canonical framework (incom- 
pressible 2D MHD with uniform resistivity) with the outside field B y $(y) 
varying on the global scale along the layer. We have first derived a system of 
rescaled equations that should be valid in the limit S ^ oo. Then, we have 
developed a 2D resistive MHD code that followed the time evolution of the 
system in order to achieve the steady state. 

We conclude that, under the assumptions of our model, the Petschek- 
like solutions are unstable and the system quickly evolves to the only stable 
steady-state solution corresponding to the Sweet-Parker reconnection layer. 
Thus, the Petschek mechanism for fast reconnection does not work in our 
model. The steady-state reconnection rate in our model problem is remark- 
ably close to the Sweet-Parker value E SP = B y ^(0)VA/\^S- 
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This main result is consistent with the results of simulations conducted 



by Biskamp|| and also those by Ugai[|12] and by Scholer|T3|. It also agrees 
with the experimental results in the MRX experiment 



Finally, even though we draw our conclusions (about Petschek-like struc- 
tures being unstable) only for this very specific model, this result is funda- 
mentally important, because this model is the canonical framework typical 
of most models of magnetic reconnection, including both Sweet-Parker and 
Petschek. This framework is the simplest possible framework for a reconnec- 
tion problem, and thus provides the necessary foundation on top of which 
one can add more complicated physical processes. Because the Sweet-Parker 
model with the classical (Spitzer) resistivity is known to be too slow to ex- 
plain the very fast time scale for the energy release in solar flares, one has to 
look for physics beyond resistive MHD with the Spitzer resistivity. The inclu- 
sion of some new physical processes into the model (for example, locally en- 
hanced anomalous resistivity is probably the most suitable candidate) would 
create a very different situation in which some Petschek-like structure with 
fast reconnection may be possible. 
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7 APPENDIX 

We give a rough physical argument for why the length V of the Petschek 
diffusion region must be of order the global length L, and correspondingly, 
why the Petschek reconnection rate must reduce to the Sweet-Parker rate. 

The equation of motion for the B x field should have two terms: a loss 
term —{va/ 'L')B x , and a source term due to the nonuniform merging of the 
incoming B y field. An estimate for the size of this source term is as follows. 

The boundary condition for the external field B y is determined globally 
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by a field scaled to the global length L. For simplicity, we take this as 

B y , cxt = B (l-y 2 /L 2 ), (Al) 

and apply it over the length of the Petschek diffusion layer, L' <C L. 

If the external field were entirely uniform, then the merging of the field 
in the diffusive region would be uniform and there would be no tendency to 
develop a B x field. If the boundary condition on B y is satisfied and V <C L, 
the tendency to develop B x by nonuniform merging would be proportional 
to L' 2 /L 2 , the nonuniformity of the external field. There are several terms 
leading to the generation of B x , all of the same order of magnitude. 

For example, consider the rate of generation of B x by the resistive diffusive 
velocity alone. Let the layer initially have only a B y field and be of uniform 
thickness 5. Then j ~ B yyCxt /5 is stronger at y = than at y — V by 
the amount Bo/5(L f2 /L 2 ), so the electric field E = rjj is correspondingly 
larger. Therefore, the rate at which lines enter at y — is larger by = 
(ri/5)(L' 2 /L 2 )B than that at y = V . 

These excess lines must turn from the y direction to the x direction and 
produce the B x field. Dividing this rate by V we get the rate of generation 
of the B x field, 

B x {generation) = (rj/5 )(L' / L 2 )B . (A2) 

One can convince oneself that the other terms, that generate the B x field 
such as the variation in the layer thickness 5, the variation in the x velocity 
of the plasma, etc., are of the same order of magnitude. Thus, the above 
rate (i]/5)(L' 2 / L 2 )(Bq/ L') should give a rough positive upper bound to the 
regeneration rate of the B x field. Therefore, the B x field cannot be stronger 
than that given by the balance equation 

dB x n V A rjL' f A ^ 

so B x = (r]/5)(L' 2 / L 2 )Bq/Va- Now, according to Petschek, the emerging 
B x field drives the shocks at a speed B x j ^Anp in the x direction, which just 
balances the incoming plasma with velocity v Tec . 
Thus, 

B x VL' 2 /A1 . 

= vm = si? {M) 

But we also have, from Petschek, t> rcc = 77/^0, so that 

L' = L, (A5) 
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and this reduces the Petschek reconnection rate 



^rec,Petschck — -^/^y \AO) 



to the Sweet-Parker rate. 
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